% The code is based on the following paper:
%   Chen, T. Y., Lin, Y. L., and Tzeng, L. Y., forthcoming.
%   Estimating probability weighting functions through option pricing bounds. Review of Asset Pricing Studies.
%
% Copyright: Tzu-Ying Chen, Yo-Lan Lin, Larry Y. Tzeng
% Date: January 30, 2024

clear

%% Setting
RootName_DIST = ['EGARCH'];

%% Load Data
FileName = ['AllTradingDate19962021_Monthly.mat'];
load(FileName, ...
     'Target_AllDate');             
clear FileName

Date_Monthly = Target_AllDate;                                
clear Target_AllDate 

%% Calculate Deviation for Each Parameters
% Transaction Cost
TC = 0.005;

% Number of Linear Segments
NUM_PL = 10;  

for d = 1:size(Date_Monthly, 1)
    Target_Date = Date_Monthly(d, 1);
    Target_TTM = Date_Monthly(d, 3);

    % Load Data 
    FileName = ['OP' num2str(Target_Date) '_TTM' num2str(Target_TTM) '.mat'];
    load(FileName, ...
         'Data', ...
         'Index_Date', 'Index_TTM', 'Index_CPFlag', 'Index_K', 'Index_OP_Bid', 'Index_OP_Ask', ...
         'Index_RF', 'Index_DY', 'Index_IS', 'Index_IS_ADJ', 'Index_IV_Bid', 'Index_IV_Ask');
    clear FileName 

    Index_Drop = find(sum(isnan(Data(:, [Index_IV_Bid Index_IV_Ask])), 2) > 0);
    Data(Index_Drop, :) = [];                                              
    AllData = Data;
    clear Index_Drop Data

    K = AllData(:, Index_K);
    S0 = AllData(:, Index_IS);
    Index_KS = size(AllData, 2) + 1;
    AllData(:, Index_KS) = K ./ S0;                                        
    AllData = sortrows(AllData, Index_KS);                                  
    clear K S0

    % Gross Return Distribution
    FileName = ['DG_GrossRET_Monthly_' num2str(Target_Date)];
    load(FileName, ...
         'RET_Monthly');
    clear FileName

    [AllRET, ~, Index_AllRET] = unique(RET_Monthly); 
    AllRET_PROP = accumarray(Index_AllRET, 1, size(AllRET), @sum, nan) / length(RET_Monthly);
    State_Monthly = AllRET';
    State_PROB = AllRET_PROP';
    clear Index_AllRET RET_Monthly AllRET AllRET_PROP

    % Slope Estimation
    K = AllData(:, Index_K);
    KS = AllData(:, Index_KS);
    OP_Bid = AllData(:, Index_OP_Bid);
    OP_Ask = AllData(:, Index_OP_Ask);
    CPFlag = AllData(:, Index_CPFlag);
    IV_Bid = AllData(:, Index_IV_Bid);
    IV_Ask = AllData(:, Index_IV_Ask);
    S0 = AllData(:, Index_IS);
    S0_ADJ = AllData(:, Index_IS_ADJ);
    TTM = AllData(:, Index_TTM) / 365;                                     
    RF = AllData(:, Index_RF);                                             
    DY = AllData(:, Index_DY);                                             

    % Expected Utility        
    OP_UB_EU = RDEU_OP_UB(State_Monthly, State_PROB, ...
                          K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                          1, 1);   

    % Rank-Dependent Expected Utility
    if sum(max(OP_UB_EU, 0) < OP_Bid') > 0
        AllTP = linspace(0, 1, 1 + NUM_PL);
        AllTP(1) = [];                                                     
        dx = 1 / NUM_PL;

        % Without Any Constraint
        TP = [];
        x0 = ones(1, NUM_PL);                                     
        A = [];                                                    
        b = [];
        Aeq = dx * ones(1, NUM_PL);  
        beq = 1;        
        LB = zeros(1, NUM_PL);           
        
        [OPT FVal ExitFlag] = fmincon(@(x) PenaltyFcn_UB(State_Monthly, State_PROB, ...
                                                         OP_Bid, K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                                         x, TP), ...
                                      x0, ...
                                      A, b, Aeq, beq, LB, [], ...
                                      @(x) NLCon_UB(State_Monthly, State_PROB, ...
                                                    OP_Bid, K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                                    x, TP), ...
                                      optimset('Algorithm', 'interior-point', ...
                                               'MaxFunEvals', 100000, ...
                                               'MaxIter', 100000, ...
                                               'TolX', 1e-10, ...
                                               'TolFun', 1e-10, ...
                                               'Display', 'Off'));             
        clear x0 A b Aeq beq LB UB
            
        OP_UB = RDEU_OP_UB_PL(State_Monthly, State_PROB, ...
                              K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                              OPT, TP);
        OP_UB = max(OP_UB, OP_UB_EU);                                      
        OP_UB = max(OP_UB, 0);                                             
        NUM_VIO = sum(OP_UB < OP_Bid');
        Record_OPT_WAC = OPT; 
        if length(TP) > 0
            Record_OPT_COND_WAC = [TP NUM_VIO ExitFlag FVal];
        else
            Record_OPT_COND_WAC = [nan NUM_VIO ExitFlag FVal];
        end
        clear TP OPT FVal ExitFlag OP_UB NUM_VIO    

        % Inverse-S Shaped
        Record_OPT = nan(NUM_PL, NUM_PL);                                  
        Record_OPT_COND = nan(NUM_PL, 4);                                  
        for k = 1:(NUM_PL - 1)
            TP = AllTP(k);
            [~, Index_TP] = min(abs(AllTP - TP));
            x0 = ones(1, NUM_PL);                                     
            EyeMatrix = eye(Index_TP);
            A_BeforeTP = - EyeMatrix + ...
                         [zeros(Index_TP, 1) EyeMatrix(:, 1:(end - 1))];
            A_BeforeTP(end, :) = [];                                       
            EyeMatrix = eye(NUM_PL - Index_TP);
            A_AfterTP = EyeMatrix - ...
                        [zeros(NUM_PL - Index_TP, 1) EyeMatrix(:, 1:(end - 1))];
            A_AfterTP(end, :) = [];                                        
            A = [[A_BeforeTP zeros(Index_TP - 1, NUM_PL - Index_TP)]; ...
                 [zeros(NUM_PL - Index_TP - 1, Index_TP) A_AfterTP]];          
            b = zeros(size(A, 1), 1);
            Aeq = [[(dx * ones(1, Index_TP)) zeros(1, NUM_PL - Index_TP)]; ...
                   [zeros(1, Index_TP) (dx * ones(1, NUM_PL - Index_TP))]];  
            beq = [TP (1 - TP)]';
            LB = zeros(1, NUM_PL);               
            
            [OPT FVal ExitFlag] = fmincon(@(x) PenaltyFcn_UB(State_Monthly, State_PROB, ...
                                                             OP_Bid, K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                                             x, TP), ...
                                          x0, ...
                                          A, b, Aeq, beq, LB, [], ...
                                          @(x) NLCon_UB(State_Monthly, State_PROB, ...
                                                        OP_Bid, K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                                        x, TP), ...
                                          optimset('Algorithm', 'interior-point', ...
                                                   'MaxFunEvals', 100000, ...
                                                   'MaxIter', 100000, ...
                                                   'TolX', 1e-10, ...
                                                   'TolFun', 1e-10, ...
                                                   'Display', 'Off'));                  
            clear EyeMatrix A_BeforeTP A_AfterTP x0 A b Aeq beq LB UB
            
            OP_UB = RDEU_OP_UB_PL(State_Monthly, State_PROB, ...
                                  K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                  OPT, TP);
            OP_UB = max(OP_UB, OP_UB_EU);                                  
            OP_UB = max(OP_UB, 0);                                         
            NUM_VIO = sum(OP_UB < OP_Bid');
            Record_OPT(k, :) = OPT; 
            Record_OPT_COND(k, :) = [TP NUM_VIO ExitFlag FVal];
            clear TP Index_TP OPT FVal ExitFlag OP_UB NUM_VIO  
        end
        Record_OPT_IS = Record_OPT;
        Record_OPT_COND_IS = Record_OPT_COND;
        clear Record_OPT Record_OPT_COND
        
        % S Shaped
        Record_OPT = nan(NUM_PL, NUM_PL);                                  
        Record_OPT_COND = nan(NUM_PL, 4);                                  
        for k = 1:(NUM_PL - 1)
            TP = AllTP(k);
            [~, Index_TP] = min(abs(AllTP - TP));
            x0 = ones(1, NUM_PL);                                     
            EyeMatrix = eye(Index_TP);
            A_BeforeTP = EyeMatrix - ...
                         [zeros(Index_TP, 1) EyeMatrix(:, 1:(end - 1))];
            A_BeforeTP(end, :) = [];                                       
            EyeMatrix = eye(NUM_PL - Index_TP);
            A_AfterTP = - EyeMatrix + ...
                        [zeros(NUM_PL - Index_TP, 1) EyeMatrix(:, 1:(end - 1))];
            A_AfterTP(end, :) = [];                                        
            A = [[A_BeforeTP zeros(Index_TP - 1, NUM_PL - Index_TP)]; ...
                 [zeros(NUM_PL - Index_TP - 1, Index_TP) A_AfterTP]];          
            b = zeros(size(A, 1), 1);
            Aeq = [[(dx * ones(1, Index_TP)) zeros(1, NUM_PL - Index_TP)]; ...
                   [zeros(1, Index_TP) (dx * ones(1, NUM_PL - Index_TP))]];  
            beq = [TP (1 - TP)]';
            LB = zeros(1, NUM_PL);               
            
            [OPT FVal ExitFlag] = fmincon(@(x) PenaltyFcn_UB(State_Monthly, State_PROB, ...
                                                             OP_Bid, K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                                             x, TP), ...
                                          x0, ...
                                          A, b, Aeq, beq, LB, [], ...
                                          @(x) NLCon_UB(State_Monthly, State_PROB, ...
                                                        OP_Bid, K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                                        x, TP), ...
                                          optimset('Algorithm', 'interior-point', ...
                                                   'MaxFunEvals', 100000, ...
                                                   'MaxIter', 100000, ...
                                                   'TolX', 1e-10, ...
                                                   'TolFun', 1e-10, ...
                                                   'Display', 'Off'));           
            clear EyeMatrix A_BeforeTP A_AfterTP x0 A b Aeq beq LB UB
            
            OP_UB = RDEU_OP_UB_PL(State_Monthly, State_PROB, ...
                                  K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                  OPT, TP);
            OP_UB = max(OP_UB, OP_UB_EU);                                  
            OP_UB = max(OP_UB, 0);                                        
            NUM_VIO = sum(OP_UB < OP_Bid');
            Record_OPT(k, :) = OPT; 
            Record_OPT_COND(k, :) = [TP NUM_VIO ExitFlag FVal];
            clear TP Index_TP OPT FVal ExitFlag OP_UB NUM_VIO 
        end
        Record_OPT_S = Record_OPT;
        Record_OPT_COND_S = Record_OPT_COND;
        clear Record_OPT Record_OPT_COND
        
        % Globally Concave
        TP = [];
        x0 = ones(1, NUM_PL);                                     
        EyeMatrix = eye(NUM_PL);
        A = - EyeMatrix + ...
            [zeros(NUM_PL, 1) EyeMatrix(:, 1:(end - 1))];
        A(end, :) = [];                                                    
        b = zeros(size(A, 1), 1);
        Aeq = dx * ones(1, NUM_PL);  
        beq = 1;        
        LB = zeros(1, NUM_PL);           
        
        [OPT FVal ExitFlag] = fmincon(@(x) PenaltyFcn_UB(State_Monthly, State_PROB, ...
                                                         OP_Bid, K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                                         x, TP), ...
                                      x0, ...
                                      A, b, Aeq, beq, LB, [], ...
                                      @(x) NLCon_UB(State_Monthly, State_PROB, ...
                                                    OP_Bid, K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                                    x, TP), ...
                                      optimset('Algorithm', 'interior-point', ...
                                               'MaxFunEvals', 100000, ...
                                               'MaxIter', 100000, ...
                                               'TolX', 1e-10, ...
                                               'TolFun', 1e-10, ...
                                               'Display', 'Off'));             
        clear EyeMatrix x0 A b Aeq beq LB UB
            
        OP_UB = RDEU_OP_UB_PL(State_Monthly, State_PROB, ...
                              K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                              OPT, TP);
        OP_UB = max(OP_UB, OP_UB_EU);                                      
        OP_UB = max(OP_UB, 0);                                             
        NUM_VIO = sum(OP_UB < OP_Bid');
        Record_OPT_CC = OPT; 
        if length(TP) > 0
            Record_OPT_COND_CC = [TP NUM_VIO ExitFlag FVal];
        else
            Record_OPT_COND_CC = [nan NUM_VIO ExitFlag FVal];
        end
        clear TP OPT FVal ExitFlag OP_UB NUM_VIO 
        
        % Globally Convex
        TP = [];
        x0 = ones(1, NUM_PL);                                     
        EyeMatrix = eye(NUM_PL);
        A = EyeMatrix - ...
            [zeros(NUM_PL, 1) EyeMatrix(:, 1:(end - 1))];                  
        A(end, :) = [];
        b = zeros(size(A, 1), 1);
        Aeq = dx * ones(1, NUM_PL);  
        beq = 1;        
        LB = zeros(1, NUM_PL);           
        
        [OPT FVal ExitFlag] = fmincon(@(x) PenaltyFcn_UB(State_Monthly, State_PROB, ...
                                                         OP_Bid, K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                                         x, TP), ...
                                      x0, ...
                                      A, b, Aeq, beq, LB, [], ...
                                      @(x) NLCon_UB(State_Monthly, State_PROB, ...
                                                    OP_Bid, K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                                    x, TP), ...
                                      optimset('Algorithm', 'interior-point', ...
                                               'MaxFunEvals', 100000, ...
                                               'MaxIter', 100000, ...
                                               'TolX', 1e-10, ...
                                               'TolFun', 1e-10, ...
                                               'Display', 'Off'));           
        clear EyeMatrix x0 A b Aeq beq LB UB
            
        OP_UB = RDEU_OP_UB_PL(State_Monthly, State_PROB, ...
                              K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                              OPT, TP);
        OP_UB = max(OP_UB, OP_UB_EU);                                      
        OP_UB = max(OP_UB, 0);                                             
        NUM_VIO = sum(OP_UB < OP_Bid');
        Record_OPT_CV = OPT; 
        if length(TP) > 0
            Record_OPT_COND_CV = [TP NUM_VIO ExitFlag FVal];
        else
            Record_OPT_COND_CV = [nan NUM_VIO ExitFlag FVal];
        end        
        clear TP OPT FVal ExitFlag OP_UB NUM_VIO
        clear AllTP dx
    else
        Record_OPT_WAC = [];
        Record_OPT_COND_WAC = [];        
        Record_OPT_IS = [];
        Record_OPT_COND_IS = [];        
        Record_OPT_S = [];
        Record_OPT_COND_S = [];        
        Record_OPT_CC = []; 
        Record_OPT_COND_CC = [];        
        Record_OPT_CV = []; 
        Record_OPT_COND_CV = [];
    end
    
    % Output
    FileName = ['Summary_DEV_' num2str(Target_Date)];
    save(FileName, ...
         'Record_OPT_WAC', 'Record_OPT_COND_WAC', ...
         'Record_OPT_IS', 'Record_OPT_COND_IS', ...
         'Record_OPT_S', 'Record_OPT_COND_S', ...   
         'Record_OPT_CC', 'Record_OPT_COND_CC', ...   
         'Record_OPT_CV', 'Record_OPT_COND_CV', ...   
         'OP_UB_EU', ...
         'AllData', ...
         'OP_Bid', 'OP_Ask', 'IV_Bid', 'IV_Ask', ...
         'K', 'KS', 'CPFlag', 'S0', 'S0_ADJ', 'TTM', 'RF', 'DY', ...
         'TC', 'NUM_PL', ...
         'State_Monthly', 'State_PROB');
    clear FileName    
    
    clear Target_Date Target_TTM AllData State_Monthly State_PROB
    clear Index_Date Index_TTM Index_CPFlag Index_K Index_OP_Bid Index_OP_Ask
    clear Index_RF Index_DY Index_IS Index_IS_ADJ Index_IV_Bid Index_IV_Ask Index_KS 
    clear K KS CPFlag OP_Bid OP_Ask IV_Bid IV_Ask S0 S0_ADJ TTM RF DY
    clear OP_LB_EU OP_UB_EU Record_OPT_WAC Record_OPT_COND_WAC
    clear Record_OPT_IS Record_OPT_COND_IS Record_OPT_S Record_OPT_COND_S
    clear Record_OPT_CC Record_OPT_COND_CC Record_OPT_CV Record_OPT_COND_CV
end
clear TC NUM_PL
clear d k